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Abstract 

Nonlinear stochastic differential equation models with unobservable vari- 
ables are now widely used in the analysis of PK /PD data. The unobservable 
variables are often estimated with extended Kalman filter (EKF), and the 
unknown pharmacokinetic parameters are usually estimated by maximum 
likelihood estimator. However, EKF is inadequate for nonlinear PK/PD 
models, and MLE is known to be biased downwards. A density-based Monte 
Carlo filter (DMF) is proposed to estimate the unobservable variables, and a 
simulation-based procedure is proposed to estimate the unknown parameters 
in this paper, where a genetic algorithm is designed to search the optimal 
values of pharmacokinetic parameters. The performances of EKF and DMF 
are compared through simulations, and it is found that the results based 
on DMF are more accurate than those given by EKF with respect to mean 
absolute error. 

Keywords: PK/PD modeling, Stochastic differential equation, Extended 
Kalman filter, Density-based Monte Carlo filter, Genetic algorithm 



* Corresponding author. 
Email addresses: hgh@cqu.edu.cn (Guanghui Huang), hust_jp_w@yahoo.com.cn 
(Jianping Wan), chenhui_tj@126.com (Hui Chen ) 



January 19, 2013 



1. Introduction 



Stochastic differential equations (SDEs) are powerful tools in pharma- 
cokinetic and pharmacodynamic (PK/PD) modeling, which can be used as 
a diagnostic tool to facilitate systematic model development H or as a re- 
alistic method to describe the variations in system 04a] • 0] documents that 
SDEs provide a more realistic description of the variability that improves 
individual simulation and predictive properties, accelerates model speed by 
simplifying inter-occasion variability, and finally changes the model into one 
that could not be falsified by the autocorrelation function. 

Pharmacokinetic parameter estimation is one of the important steps in 
PK/PD data analysis. Maximum likelihood estimation (MLE) based on the 
extended Kalman filter (EKF)is usually applied to estimate the parameters 
in SDE models, such as [lrH,0-S] and the references therein. On the other 
hand, Kalman filter is designed to estimate the state variable involved in 
a linear model, and EKF is the linearized version of Kalman filter for the 
models with nonlinear characteristics. [3] documents that the failure to 
produce Gaussian residuals with EKF may indicate which is inadequate for 
nonlinear modeling in PK/PD data analysis, possibly motivating the pursuit 
of higher order filters or other estimation methods. 

[9| argues that even if the higher-order nonlinear filters deduced from 
Kalman filter give us less biased filtering estimates than the EKF, the fil- 
tering estimates obtained from the higher-order nonlinear filters are still 
biased because the nonlinear functions are approximated ignoring the other 
higher-order terms. And the other filters, such as the density-based Monte 
Carlo filters (DMF), might be less biased than the EKF, as the unobserv- 
able variable can be generated from the nonlinear functions directly without 
approximations. The purpose of this paper is to compare the performances 
of EKF and DMF under a nonlinear model, and to develop more efficient 
algorithms for PK/PD parameter estimation. 

As MLE is often inefficient and biased for finite sample (lol . 11]. a 
simulation-based procedure is proposed to estimate the unknown param- 
eters in this paper, where a genetic algorithm is designed to search the 
optimal values of parameters. A one-compartment pharmacokinetic model 
with nonlinear absorption and first order elimination is used to compare 
the performances of EKF and DMF through simulated investigations. It is 
found that the results based on DMF is more accurate than those based on 
EKF with respect to mean absolute error. 

The remainder of the paper is constructed as follows: Section 2 intro- 
duces the model to be investigated, Section 3 gives the EKF algorithm for 



this model, and Section 4 demonstrates the proposed DMF algorithm. The 
estimates of the unobservable variables by EKF and DMF are compared in 
Section 5. The criterion of estimation for the unknown parameters and the 
genetic algorithm of optimization are given in Section 6. The conclusions 
and discussions are given in Section 7. 



2. Stochastic nonlinear model 

A one-compartment model with nonlinear absorption and first order 
elimination is considered in this paper, which is used to describe the PK 
of a drug following an oral dose by where pharmacokinetic parameters 
are estimated with MLE based on EKF. 

Let Q(t) (mg) be the amount of drug in the GI tract at time t, t € (0, T] 
(min), which is an unobservable variable. C{t) (mg/1) is the concentration 
of drug in plasma, which is an observable variable. A system of stochastic 
differential equations are used to describe the processes of absorption and 
elimination of the drug, i.e. 

Vmax Q(t) 

K m + Q(ty 
v max Q(t) c L c(ty 



d Q(t) = - T ™„Z dt + a q dB(t), (1) 



dC(t) 



dt + a c dW{t), (2) 



(K m + Q{t)) V V 

where V max (mg/min) is the maximum reaction rate, K m (mg) is the Michaelis 
constant, a q and a c are the diffusion parameters of Q(t) and C{t) respec- 
tively, Cl (1/min) is the rate of elimination, and V (1) is the apparent vol- 
ume of distribution. The two stochastic processes B(t) and W(t), t £ [0,T], 
are two independent standard Wiener processes starting from zero. Let 
d = (V-max, K m ,V,CL, a 2 q ,al), which is a vector of six elements, and 9 £ 
0Cl 6 , where is the parameter space. The purpose of this paper is to 
estimate 6 from the limited observations of C(t), {c^ , cj 2 , • • • ,Q n }, where 
{ti,t2, ■ ■ ■ ,t n } are the time points of observations, and n is the number of 
observations. 

SDEs (pQ) and ([2]) are nonlinear equations, it is difficult to find explicit 
solution for such SDEs. In order to simulate Q(t) and C(t) at discrete time 
points, ([I]) and ([2]) are approximated with discrete differences in ltd type, 



i.e. 



Qk — Qk-i — ™ ax ^ 1 (tk — tk-i) + cr q (B k — -Bfc-i) ) (3) 
-ft-m + l°ik-l 



Ck — Ck-i + 



Vmax Qk-l Cl Cfc-l 



(iT m + Q fc _i) V V 
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where Q k = Qt k , C k = C tk , B k = B tk and W k = W tk , k = 1,2, ••• ,n. 
B(t) and W(t) are two independent standard Wiener processes, such that 
the increments B k — B k -i and W k — W k -\ are independent and identically 
distributed, where B k — B k _\ ~ N(0, ^Jt k — tt-i). 



3. Extended Kalman filter 



Let c k = ct k , which is the observation at time t k . And the information 
set at time t s is Y s = {ci, C2, • • • , c s }, where s G {1, 2, • • • , n}. Let 



(5) 



which is the conditional expectation of Q k given the information set at time 
t s . When s < k, (J5]) is called the prediction of Q k ; when s = k, ([5} is called 
the filtering of Q k ; and © is called the smoothing of Q k when s > k. 

Kalman filter is particularly powerful and useful for the linear models 
which include unobservable components. Applying the linearized nonlinear 
functions to the Kalman filter, the resulted algorithm is called the extended 
Kalman filer (EKF). Approximate the two nonlinear functions in ([3]) and 
([4]) with first order Taylor series expansion, we have the following EKF 
algorithm 



sec 



Appendix A): 



Qk\k-1 



Qfe-iiit-] 



Vmax Qk-l\k-l 
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-i 
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(6) 
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(t k — t k -i) , (8) 

(9) 
(10) 

(11) 
(12) 
(13) 



where 



J k\k-1 



Fk\k-i 



V K 

v max - £i r? 



{K m + > 



Jk-l\k-l) 

V K 

v max lv n 



V 



(tk — tk-l) 



(Km + Q 



(tk — tk-l) ■ 



(14) 
(15) 



fc-ilfc-ij 
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Set Q |o = Qoj S |o = 0, Co = 0, the unobservable variable Qk can be 
estimated by EKF algorithm in a recursive manner. 



4. Density-based Monte Carlo filter 

Density-based Monte Carlo filter is an alternative solution to nonlinear 
filtering problems, and the resulted algorithm is easy and convenient to 
compute the filtering estimate Qk\k 01 • The filtering estimation based on 
Monte Carlo technique is given by 



N 



Qk\k — Qj,k ^~>j,k, 
3=1 



(16) 



where Qj k is the simulated value of the unobservable variable at time tk in 
the jth path, which is generated from equation ([3]) directly, and ./V is the 
number of simulated paths, ojjk is the weight of the jth path at time tk, 
which satisfies ujjfi = 1/N and 



N 
3=1 



1. 



(17) 



bjjk is calculated with a recursive formula 



_ P (Cfc|Cfc-l, Qj,k-l) ■ Uj,k-1 
Z^i=l P { c k\Ck-l, Hj,k-l) • w i,fc-l 



(18) 



where p (ck\ck-i,Qj : k-i) is the conditional density function of Ck given by 
(SD, i.e. 



where 



P(Ck\Ck-l,Qj,k-l) 



2vr Ok 



cxp 



(c fc - m k y 



(19) 
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0"c (tk — tk-l) , 



Cfc-1 + 



Vmax Qj,k—1 Cl Cfc— 1 



(K m + Q jtk _ l ) y 



(tk — tk-x) ■ 



(20) 
(21) 



Details of DMF can be found in Appendix B 
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5. Estimates of drug in GI tract 



Qk is the amount of drug in the GI tract at time t^, which is an unob- 
servable variable when tk £ (0, T]. The concentrations of drug in plasma can 
be observed at different time points, and Qk can be estimated from those 
observations by EKF and DMF respectively. In order to compare the per- 
formances of EKF and DMF, a simulated investigation is designed in this 
paper. 



5 mg, c 



= mg/1, C L = 0.05 1/min, V 
J2 — n nnm ~2 



5 1, V„ 



1 



Set Q 

mg/min, K m = 15 mg, a l q = 0.0002, a l c = 0.00003, and t = 5, 10, 15, 20 
25, 30, 40, 50, 60, 90, 120, 150, 180, 230, 290, 340, and 390 min. There 
are 17 observation time points, and the corresponding amounts of drug Qk 
and concentrations in plasma Ck are generated from equation ([3]) and (j4]) 
respectively. The filtering estimate of Qk is Qk\ki which is calculated from 
those observed values of concentrations by EKF and DMF respectively. A 
plot of the observed C(t) and the estimated Qk\k versus time is given in 
Figure [TJ 



Figure 1: The observed concentrations in plasma and the estimated amounts of drug in 
GI tract. 
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Table 1: MAEs of estimates for drug in GI tract by DMF and EKF within 200 simulations. 



Quantiles 0.0500 0.3000 0.5000 0.6000 0.7000 0.8000 0.9000 0.9500 

DMF 0.0233 0.0335 0.0399 0.0418 0.0448 0.0487 0.0546 0.0591 

EKF 0.0600 0.0738 0.0797 0.0825 0.0863 0.0908 0.0990 0.1021 

RD 1.5747 1.2042 0.9997 0.9735 0.9246 0.8647 0.8117 0.7288 



The mean absolute error (MAE) is denned as 

1 n 

MAE = -Y,\Qk-Qk\k\, (22) 

which is used to measure the accuracy of estimates, where n is the number 
of observations. 

The simulated investigation is repeated 200 times, and the MAE of each 
simulation is calculated. A quantile analysis is applied to those observed 
MAEs, and the results are reported in Table [H where DMF and EKF in- 
dicate the results are given by DMF and EKF respectively. And RD is the 
relative difference between the values of DMF and EKF, i.e. 

MAE of EKF - MAE of DMF 
RD ~ MAE of DMF • (23) 

It is found that the 0.95 quantile of MAEs given by DMF is smaller than the 
0.05 quantile given by EKF. It can be concluded that the errors of estimates 
given by DMF is much smaller than their counterparts given by EKF. This 
result can be regarded as another evidence to support the argument in [jj, 
where EKF is found to be inadequate for nonlinear modeling in PK/PD 
data analysis. 



6. Estimation of parameters 

6.1. Criterion of estimation 

MLE based on the extended Kalman filter (EKF) is used to estimate the 
parameters in SDE models by several authors, such as and the 

references therein. On the other hand, MLE is often inefficient and biased 
for finite sample. As the sample size is limited in this paper, an alternative 
criterion of estimation is adopted to estimate the unknown parameters. 
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For a particular parameter 9, the filtering estimate of Qk is denoted 
as Qfeifc(#), which is a function of 9. Substitute Qk~i with Qk-i\k—l(@) i n 
equation ([!]), and simulate M observations of Ck from this equation, denoted 
as c^fc, c 2) fc, • • • ,0^. Let 

M 

p(c fc , 0) = > y | c k - Cj )fe |, (24) 

3=1 

where is the observed value of Ck at time tk- Let 

ri 

L(0) = X>(c*,0), (25) 
fe=l 

which is the loss function to be used in the following sections. The parameter 
9 which satisfies 

9 = arg min L(9) (26) 

is used as the estimator of 9* which generates those observed data. 9 is a 
simulation-based quasi-robust estimator, which is insensitive to departures 
from underlying assumptions. 

6.2. Optimization procedure 

The objective function (125|) is a nonlinear function, where Cj k is simu- 
lated from equation @ based on the filtering estimate of Qk, such that 9 
can not be computed explicitly. In order to solve the nonlinear optimization 
problems in PK/PD data analysis, a quasi-Newton method based on BFGS 
updating formula is adopted by several authors, such as [H,[l3] an d the refer- 
ences therein, where the gradient of the objective function is approximated 
by a set of finite difference derivatives. This algorithm can be shown to 
converge to a possible local minimum [l~2l ]. 

In order to avoid the attraction of local minimum, a genetic algorithm 
is proposed in this paper. The steps of this algorithm are as follows: 

1. Start. Generate random population of S parameters 9i £ , where 
i = l,2, -- ,S. 

2. Fitness. Evaluate the fitness of each parameter in the population with 
L(9), where the smaller L(9), the better fitness. 

3. New population. Create a new population by repeating following steps 
until the new population is complete. 
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(a) Selection. [aS pop ] parameters are selected from the population 
according to their fitness, [x] indicates the largest integer which 
is less than x, S pop is the size of the current population, and a is 
the proportion of parameters selected to be new population. 

(b) Crossover. Denote the selected parameters 61,62,- ■■ ,9c, which 
are sorted in increasing order according to their values of L{6). 
The crossover probability for the ith parameter is chosen to be 

Pi = (27) 

where i = 1, 2, ■ ■ ■ , C — 1, and Lj = L(0j). 61,62, - ■ ■ , 9c-i are 
randomly selected to cross over according to those probabilities. 
Suppose 6i and 6j are selected to cross over in the rth run, where 
Li < Lj, r = 1, 2, • • • , R, and 0j = (9 it i,9 ij2 , ■ ■ ■ , 6 ifi ). Two new 
parameters are generated in the following way, 

Q™ ew = 0i + (9i - 6j) x weight x temperature, (28) 
Qnew _ q,_ ^. _q.^ x we igfo x temperature, (29) 

weight = — — (30) 
Lj + Lj 



where temperature is a parameter to control the speed of con- 
vergence. If the jth element 6^ w of the new parameter 6™ ew = 
(6™!™, 9^2", • • • , 9™f") is lager than the upper bound 9j, or smaller 
than the lower bound 6j , then set 

qnew _ 36jj + 6j 39jj + 9j 

9 r ,j 1 . or 6 rj . (31) 

Which is also true for 9™ ew ■ 

(c) Mutation. Another S parameters are generated randomly from 
the set 6, denoted as 9f, 9f,- ■ ■ , 9f. 

(d) Accepting. Evaluate the values of loss function at 6 1 , 62, ■ ■ ■ , 6c, 

anew anew anew anew anew anew aM qM qM „„j 

select the front [a(C+2R+M)] parameters as the new population. 

4. Replace. Use new generated population for a further run of algorithm. 

5. Test. Denote the populations in the ith and (i + l)th generations as 
{9\, 9\, ■ ■ ■ , 9* ni } and ■ • ■ , ^+ + \ } . Set a series of proba- 
bilities < ct\ < 02 < • • • < a a < 1 and find those corresponding 
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quantiles of these two populations, denoted as {0 l ai ,6 l a2 , • • • ,6 l a } and 

{e i +\e i +\---,9 i + 1 }.Let 

et , - ir 1 . . 



the end condition is chosen to be EC < 0.00001, or the number of 
loops is beyond 100. If the end condition is satisfied, stop and return 
the best solution in current population. The estimator 6 is the mean 
of the last population. 
6. Loop. Go to step 2. 

Let 0\ and be the first parameters in the ith and (i + l)th genera- 
tions, which must satisfy 

o < l (e[ +1 ) < l (e[) < l (el) , (33) 

which ensures the convergence of algorithm. The mutation step reduces the 
risk of attraction of local minima. The proposed genetic algorithm needs 
not to approximate the gradients with a set of finite difference derivatives, 
and the burden of programming is much less than the quasi-Newton method 
based on the BFGS updating formula. 



1 6 a 

7=1 k=i 



6.3. Estimates comparison 

Set Q = 5 mg, c = mg/1, C L = 0.05 1/min, 7 = 51, V max = 
1 mg/min, K m = 15 mg, a\ = 0.0002, a\ = 0.00003, and t = 5,10, 
15,20,25,30,40, 50,60,90,120,150, 180,230,290,340, and 390 min. There 
are 17 time points, and those observations are generated from ([3]) and ([!]). 
Parameter 9 = (V max , K m , V, Cl, o~g, a^.) will be estimated from the limited 
observations by the proposed estimators based on EKF and DMF respec- 
tively. 

The simulated experiment is repeated 200 times in this paper, and the 
vector of parameters is estimated in each simulation, where temperature = 
0.75, the proportion a used to select new populations is 0.05, and the quan- 
tiles used to construct the end condition are 0.2, 0.4, 0.5, 0.6, and 0.8. A 
quantile analysis is applied to those estimated parameters, and the results 
are reported in Table [2j There are 11 quantiles reported in this table, in- 
cluding 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9 and 0.95. The real value 
of each parameter is given in the line Real. And entries in the first row of 
each quantile are the estimated parameters given by the algorithm based on 
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DMF, and the second row are those estimates given by the method based 
on EKF respectively. 

The mean absolute error of estimated parameters is defined as 

1 N 6 

MAEP = ^££l^-*H ( 34 ) 
i=i j=i 

where 6ij is the jth element of the ith estimated parameter in N simula- 
tions, and 9* is the jth element of the real parameter which generates the 
observed data. The MAEP of EKF is 0.7070 among 200 simulations, and 
its counterpart of DMF is 0.5929, where the later is 16.13% less than the 
former, which indicates that the estimator based on DMF is better than the 
one based on EKF with respect to MAEP. 

The algorithms proposed in this paper are programmed with Matlab 
R2009, which run on a personal computer with an Intel(R) Core(TM)2 Duo 
CPU E7500, whose main frequency is double 2.93 GHz. 

7. Conclusions and discussions 

A density-based Monte Carlo filter is proposed to estimate the unob- 
servable variables in a nonlinear stochastic differential equation model, and 
a simulation-based quasi-robust estimator is proposed to estimate the un- 
known pharmacokinetic parameters in this model. A genetic algorithm is 
proposed to solve the optimization problem in the estimation procedure. 
The performance of the proposed filter is compared with the extended 
Kalman filter, and it is found that DMF is more efficient than EKF in 
the simulation investigations. 

Further research possibilities are mainly in three directions. First of 
all, other nonlinear filters can be applied in the analysis of PK/PD data. 
Several nonlinear filters are used to estimate the unobservable variables in 
state-space models, including the Gaussian sum filter, the numerical inte- 
gration filter, the importance sampling filter, the rejection sampling filter, 
and the density-based Monte Carlo filter. It should be possible to determine 
the optimal filtering algorithm for a particular PK/PD model. The second 
direction concerns the estimation criterion that can be used in the analysis 
of PK/PD data. It is often stated that MLEs are biased for finite sample, 
while the sample size in PK/PD data analysis is often limited. The third 
direction is the algorithm to be used in the procedure of optimization. The 
efficiency of optimizing algorithms should be taken into account in PK/PD 
data analysis. 
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Appendix A. Algorithm of EKF 



Suppose y k is the value of observable variable at time t k , and a k is the 
unobservable state variable at time t k , which satisfy 



y k = h k (a k ,e k ), 
a k = g k (a k -i,r] k ), 



(A.l) 
(A.2) 



where e k and i] k are two independent disturbances at time t k . h k (-) and 

g k (-) are two nonlinear functions, which can be approximated with first 
order Taylor series expansions 

y k pa h k \ k _x + Z k \ k _ x (a k - a, k \ k -\) + S k \ k _ie k , (A. 3) 

«fc ~ 9k\k-l + Tk\k-i( a k-i — ajfc-l|fc-l) + Rk\k-iVk, (A. 4) 



where 



l k\k-l 



J k\k-i 



5 fc|fc-i 



9k\k-l 



l k\k-l 



^fe(afe|fc-i.0), 
dh k (a k ,e k 



da k 
dh k (a k ,e k ) 



{(Xk,£k)=(<Xk\k-lfi) 



d £ k K,e t )=K|t-i-0)- 

3fc(«fe-i|fc-i,0), 
dg k (a k -i,r)k) 



R 



k\k-i 



da k -i 
dg k {a k -i,r)k) 



dm 



(ak-i,r]k)=( a k-i\k-ifi), 
(afc_i,7/fc)=(a fc _i| fc _i,0). 
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EKF is given by the following algorithm: 

a k\k-i = 9k\k-i> (A-5) 

= r fc | fc _ 1 S fc _ 1 | A ._ 1 T^| A ,_ 1 + Rk\k-iQkR'k\k-v (A-6) 

Vk\k-i = h k \ k _i, (A. 7) 

Fk\k-i = ^fc|fe-i^fc|fe-i^fc|ft_i + Sk\k-\.HtSk\k-\i (A. 8) 

= ^fc|fc-iSjfc|fc_i, (A. 9) 

= -^fc|jfc-i-^fc[fc-i> (A.10) 

s fc|fc = s fc|fc-i - K kF k \ k ^iK' k , (A. 11) 

Ofc|fc = «fc|fe-i + K k (y k - y fe | fc _i) , (A.12) 

where £ |o = 0, a | = Qo, and Q k = H k = t k -tp-i for = 1,2, ••• , n in 
this paper. The details of EKF can be found in 



Appendix B. Algorithm of DMF 

Denote the collection of state-vector as 

M = {?o,5l,--- >?t}> (B- 1 ) 

where % is the value of unobservable variable at time t&, A; = 1, 2, ••• , n. 
The joint density function of (Yt, At) is 

p(Y t ,A t )=p(A t )p(Y t \A t ), (B.2) 

where p (At) and p (Y t |A t ) are 



p(A t ) = p(qo)J[p(q s \q s -l), (B.3) 

8=1 

t 

p(y t |A t ) = n^(c s |g s _!) , (B.4) 

s=l 

where p (g s |g s _i) and p (c s |g s _i) are obtained from ([3]) and (HJ respectively. 
The filtering density function is given by 

J p(^ s )p(y s |A s )dA s 
such that the filtering estimate of the state variable is given by 
n _ F rn Jqtp(Yt\At)p(At)dA t 

Qt\t - &[Qt\Yt\ - , 777T-T\ (A\AA ' 

J p (Yt\A t ) p (At) dA t 
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Generating random draws of At from p (At), the filtering estimate based on 
the Monte Carlo technique is give by 



^T,liQi,tP(Y t \A ht 



Q 



t\t 1 v^TV 



i?rup(Yt\A,t) 

E£=i Qi,t 111=1 p(ys\u,s-i) 



(B.7) 



Eili UUiPivslAs-i) 

where A^t is the collection of random draws for the ith generated path, i.e. 

A,* = {Qi,o, Qi,i, • • • j Qi,t} • (B.8) 



Denote 



w . IIa=lP(j/fllQj,«-l) 

J ' E5Lin*=iP(i/«i0i,«-i) 



then we have 



where 



P (2/fc|Qj,fe-i) ^i.fc-i 



Ej=i p (yfc|<9j,fc-i) Wj,fc_i 



AT 

= 1, 

and Wj )0 = l/N. The DMF of Q k is given by 



iV 

Qk\k = X] Qj,k u j,k- (B-9) 



The details of DMF can be found in 
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Table 2: Quantile analysis of parameters estimated by DMF and EKF within 200 simula- 
tions. 
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